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Abstract 

Phase-plane stability analysis of a dynamical system describing the Universe as 
a two-fraction fluid containing baryonic dust and real virial gas quintessence is 
presented. Existence of a stable periodic solution experiencing inflationary periods 
is shown. A van der Waals quintessence model is revisited and cyclic Universe 
solution again found. 



1 Introduction 



There has been an extensive search for new forms of matter, energy and possible new 
physics beyond the four fundamental forces described by the Standard Model and Gen- 
eral Relativity to conform with the present day observational picture in cosmology and, 
in particular, to account for the cosmic acceleration [1]. A fifth clement, quinta essen- 
tia, in addition to baryons, leptons, radiation, and cold dark matter, has been sought, 
typically, in the form of cosmological constant, scalar fields, hybrid models involving 
both fluid and scalar field, modification of gravity, quantum gravity models, etc — see, 
for example, [2], [3], [4], [5], [6], [7], [8], [9], [10] and the references, therein. 
A classical physics model with a two-phase fluid, described by a van der Waals equation 
of state, and with no need of introduction of scalar fields, has also been studied [11], [12], 
[13] as quintessence. In this model, negative absolute temperature is needed to explain 
accelerated expansion if the present day energy density of the dark matter (described 
as the van der Waals fluid) is much smaller than the van der Waals critical density. 
This paper presents another classical physics model describing the Universe as an infi- 
nite, flat, two-fraction mixture of baryonic dust and a real gas with equation of state 
derived from the virial expansion. Dynamical phase-plane analysis (with Hubble's pa- 
rameter H and density of the real gas p as dynamical variables) of such model shows 
that there is initial data leading to a cyclic Universe solution that goes through an 
inflationary phase in each cycle, together with open trajectories in the phase plane that 
may or may not pass through regions characterized by inflation. The temperature in 
this model is also treated as a parameter. However, this parameter that does not have 
to take negative values. It is shown that, as the Universe cools down, the inflationary 
region on the phase plane decreases and eventually disappears in the limit T — >■ 0. The 
cosmological model presented here also does not exhibit an endless sequence of cycles of 
expansion and contraction. The trajectory of the Universe on the phase place is, in first 
approximation, an ellipse and the frequency of oscillations decreases, while the ratio of 
of its axes decreases as the Universe is cooling. Eventually periodicity is lost as well. 
It is essential to point out that the particle interaction of this model has been described 
as two-particle interaction between indestructible "hard spheres" surrounded by finite 
potential wells — a regularized form of a classical real gas interaction which is slightly 
attractive at long distances and strongly repulsive at short ones. The van der Waals 
model is also revisited in the light of seeking similar cyclic solution (which is found for 
negative temperatures) . A stable- node equilibrium solution for positive temperatures is 
also present. 



2 Cosmological Setup 

The setting for the analysis of the two-fraction Universe is the Friedmann-Robertson- 
Walker-Lemaitre (FRWL) cosmology [14] with metric: 



ds'^ = g^^dx^dx" = c^dt^ - a^{t) 



1 — kr"^ 



+ r^{de'^ + sm'^ 9 dcf)'^) 



(1) 
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where a{t) is the scale factor of the Universe and k is the spatial curvature parameter. 
Einstein's equations are: 

G^,^ + Ag^^ = kT^^, (2) 

where k = SttG/c^ (in Planck units, k = 1) and the matter energy-momentum tensor 
is given by: 

T^u = {pc^+py^S^^d^,-pg^,, (3) 

with p and p — the cumulative density and pressure for both fractions. 
Friedmann equations are [15]: 



AttG , _ 3p , 1.9 / , X 

a = —{p+^)a+-Ac^a (4) 



b? = ^— pa^ + -Ac^a^ - (?k. (5) 

Only the case of flat spatial three-sections (k = 0) and without cosmological constant 
(i.e. A = 0) will be of interest. 

The Hubble parameter H is H = a/a and this will be one of the dynamical variables. 
In terms of in Planck units, with = and A = 0, Priedmann equations (4)-(5) 
become: 

H"" = \{Pb + p). (6) 

H = -liPb + P + p). (7) 
The continuity equation for the real gas 

p+ip + p)^=0 (8) 

becomes 

p + 3H{p + p) = 0, (9) 

while the continuity equation for the pressureless baryonic dust is: 

Pb + 3Hpb = 0. (10) 

Differentiating (6) with respect to time and substituting into it H from (7), p from (9) 
and ph from (10), leads to an identity. Thus, equation (6) is just an integral of equations 
(7), (9), and (10). As it can be obtained from equations (6), (7) and (9), equation (10) 
will be droppped. 

Expressing the baryonic energy density pb from equation (6) and substituting it into 
equation (7) gives the dynamical equation 

A = -Ih^-\p. (11) 

The other dynamical equation is (9): 

p = -3H{p + p), (12) 
with p being the second dynamical variable. 
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3 Real Virial Gas Quintessence 



The virial expansion relates the pressure p oi a real gas to the particle number A^, the 
temperature T and the volume V of the gas [16]: 



NksT 



N /N\2 1 

l + -F(r) + (-) G{T) + ...]. (13) 



The "departure" from ideal gas (p = NksT/V) is due to the interactions between the 
gas particles which, also, no longer have zero volumes. These are included in the above 
expression via the correction term F(T), describing the two-particle interactions; the 
correction term G(T), for the three-particle interactions, and so forth. 
The correction term F(T) is [16]: 

oo 

F{T) = 27r [\l- e~^l r'^dr, (14) 



where V{r) is the two-particle interaction potential. 

For the purposes of this work, considered here is a classical real gas model of particle 
interaction that never involves more than two particles (ignoring collisions of 3 particles 
or more amounts to disregarding the term G{T) and the subsequent terms in the virial 
expansion) and which is slightly attractive at long distances and strongly repulsive at 
short ones (the dashed curve on Figure 1). Such interaction is often presented in a 
regularized form as an ensemble of identical impenetrable "hard spheres" of radius a 
(namely, two real gas particles not being able to get closer to one another at distance 
between their centres smaller than 2a), surrounded by square potential wells of width 
ad {d > 1, d = 1 corresponds to absence of well) and depth — e (where e > 0). The 
regularized potential (the solid curve on Figure 1) of the two-particle interaction is: 

oo , < r < a, 
V{r) = { -e, a<r<ad, (15) 
, r > ad. 

It is important to point out, that representing the interactions as a collection of "hard 
spheres" only (d = 1) is not sufficient for the analysis that follows. Namely, the attrac- 
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tive part (the square well) is essential and thus d > 1. 




Figure 1: The interaction potential (dashed curve) and its regularized form (solid curve). 

With the expansion of the exponent up to and including the term proportional to 
l/{kBT), the correction F{T) becomes: 

F{T) = B-^ (16) 

and the virial expansion (13) becomes the "van der Waals" equation [16]: 

In the limit NB/V — )■ 0, this equation reduces to the usual van der Waals equation 
(which was presented before the real gas virial expansion form (13) and is not "derivable" 
from it) [16]: 



N\2 



N ^ NksT 



Phase-plane stability analysis of a dynamical system describing the Universe as a two- 
fraction fluid containing baryonic dust with energy density pf, and pressure pb = 0, and 
van der Waals gas quintessence will be discussed later in this paper. Firstly, this analysis 
will be done for a quintessence of real gas whose two-particle interaction potential is 
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described by (15). In Planck units (fc^ = 1), the correction term F{T) for such gas is 
given by: 



F{T) = 27r 



ad „ 



J r^dr + J r^(l — eT^dr 



[1 + (1 - eT)(c?3 _ 1)] = -az{T), (19) 



where z{T) = {ei" - l){d^ - 1) - 1 and a = {2/3)'Ka^. 

For gases containing spherical molecules only, with a square well interaction potential 
(15) between them, the next virial coefficient G{T) is also known [17] — it is a term 
cubic in er — 1. However, the inclusion of such asymptotically smaller term does not 
significantly alter the physics of the model presented here and for that purpose three- 
particle interactions or interactions involving collisions of four or more particles will be 
disregarded. 

In terms of the gas density, the real gas equation of state is: 

p = Tp[l-apz{T)]. (20) 

In the analysis, a, e, and d will be the parameters of the model. The temperature 
T will also be treated as a parameter, but one that will be allowed to vary, and the 
different values of this parameter, would single out different stages of the "evolution" of 
the trajectories in the phase-plane of the dynamical system. Upon substitution of the 
equation of state (20), the dynamical system becomes: 

p = -SHp[l + T-apTziT)] = fi{p,H), (21) 
H = -^H^-^Tp[l-apz{T)] = h{p,H). (22) 

The first observation is that the origin of the phase-plane is an equilibrium point — at 
this point, p = = H. Equilibrium points are standardly denoted with stars and the 
equilibrium point at the origin is therefore: (p^ = 0, = 0). 

The next observation is that for equilibrium points to occur away from the origin of the 
phase-plane, the pressure p must not be positive. This is visible from equation (11) — 
at an equilibrium point, if = 0, thus, if H* / at that point, then (11) is consistent 
only when p < {H has to be real). As a is positive (it is one half of the volume of the 
"hard spheres" ) and the density p is also positive, the pressure p is negative when: 

1 - apz{T) < 0. (23) 

Therefore, equilibrium points, for which H* ^ 0, can occur only for values of the density 
p greater than 

^'"^ = az{f) = «[(e^-l)(d3-l)-l]- ^^^^ 

Clearly, z{T) = (er — l)(d^ — 1) — 1 must be positive in order pmin to be positive. This 
puts an upper limit on the temperature below which equilibrium points, different from 
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the origin of the phase-plane, exist: 



ln#-ln(d3-l)- ^^^^ 

This temperature is known from the theory of molecular gases as Boyle temperature 
Tb [17]. It is called Boyle temperature because at around this temperature, Boyle's law 
holds best [17]. This is the temperature at which the first virial coefficient F{T) — or, 
for our purposes, z(T) — changes sign. In our analysis, according to the inequality (23), 
the temperature of the virial gas must stay at all times bellow the Boyle temperature. 
This is the reality condition for the Hubble parameter H, if equilibrium points are to 
exist away from the origin of the phase plane. 

For example, the Boyle temperature is 123 K for Ne, 410 K for A, 594 K for Kr, 772 K 
for Xe [17]. 

Relative to the Boyle temperature Tb, one can define "van der Waals molecular volume" 
vb diS follows: 



VB 



dT 



T=Tb ^"^^^ 



— as it is equal to the term interpreted as molecular volume in the van der Waals 
equation of state [17]. The relationship between a and vb is [17]: 

— = (f\\nd^-\n{d^-l)]. (27) 

a 

The van der Waals volume vb is 35.5 for Ne, 67.4 for A, 78.5 A^ for Kr, 114 A^ 
for Xe [17]. 

For the rare gases, the value of d is about 2.0 [17]. 

It will be shown further in the analysis that a cyclic Universe scenario with an infla- 
tionary phase emerges within this set-up. Sufficient expanding phase can be provided 
by choosing parameters of the model in such way that the Boyle temperature is suffi- 
ciently high. For example, Boyle temperature of the order of 10^° K can be achieved by 
increasing the range of the interaction by a factor of 1000, while keeping e fixed. Also, 
increasing e increases linearly the Boyle temperature. 

For a fixed Hubble parameter H [i.e. = in (22)], the temperature T increases mono- 
tonically with the density p and tends to the Boyle temperature as p grows to infinity. 
The domain of the function T{p) is from p^^^ [see (24)] to oo. There is a minimum of 
the function outside this domain [H is not a real number there): at point /9min/2 and 
the value of the temperature at this point is the solution of the transcendental equation 
^H'^ = T/[Aaz{T)]. 

Returning to the phase-plane analysis, it can be easily seen that there is another equi- 
librium point for which p2 = Pmin , = 0. There are two more equilibrium points. 
From p = 0, it is straightforward to find 
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Substituting this into (22) and requesting H = 0, leads to: 



H*2,3 = = ± ^JIpU1 + ^). (29) 

The last two equilibrium points arc therefore (p*^ , and {p*^ , HS^). 

Next, the dynamical system is linearized near an equilibrium point (p*,H*): 

p ^ h{p,H) = {^y{p-p*) + {^yiH-H*) + ..., (30) 

H ^ Mp,H) = {^yip-p*) + (^y(H-H*) + ..., (31) 

where the stars on the derivatives indicate that they are taken at the equilibrium point 
{p*,H*). In matrix form this can be written as: 

j^X{t) = L{p*,H*).X{t), (32) 

where: 

/ P{t)-P* 

Xit) = 

\ H{t) - H* 
The stability matrix L(p, H) is: 



dfi 


dfi 


dp 


dH 


a/2 


9/2 


dp 


dH 



and L{p,H)=\ \. (33) 



r(nm- -^^^^ + ^) + ^HaTpz{T) -3p[l + T - aTpz{T)] \ 

To determine the type of an equilibrium point, the eigenvalues of the stability matrix 
should be determined at that point. 

Firstly, at the origin (p| = , Hf = 0) , the characteristic equation of the stability 
matrix L{0,0) is = 0. Thus, A = is a double eigenvalue and to analyze this 
equilibrium point and see if the origin is stable under perturbations, a small infinitesimal 
perutrbation parameter S > will be introduced in the second equation of the dynamical 
system (21)-(22): 

p = -3Hp[l + T -Tapz{T)], (35) 
H = -lH' + ^-d'-lTp[l-apz{T)]. (36) 

The introduction of the (5-term leads to a bifurcation of the original equilibrium point 
{pI = ,Hl = 0) into a "dipole" ot two equilibrium points P'^ and P~ with coordinates 
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on the phase plane {pi 
21. 



0,H[* 



+5) and (pj 



, H'{* 



-5), respectively (Figure 




Figure 2: The phase-plane diagram of the perturbed system. Points R and S are saddles, Q — 
centre, point — stable node, point P~ — unstable node. 

At the upper point of the "dipole", both eigenvalues are negative: Ai = —35(1 + r) < 
and A2 = —36 < 0. This point therefore corresponds to a stable node. At the lower 
point of the "dipole", both eigenvalues change only their signs. This is an unstable 
node as both eigenvalues are positive. Upon taking the limit (5 — )• 0, the stable and the 
unstable node coalesce and thus the origin [pi = , HI = 0) will attract trajectories in 
the phase-plane, as shown on Figure 2. 

Next, at the equilibrium point R with coordinates (^3 , //|) (determined in (28)-(29) 
above), the eigenvalues of L{p^ ,-^2) ^-^^^ 



Ai = ^3pl (1 + T) > 0, 
A2 = --fiff^ < 0. 



(37) 
(38) 



Therefore, as the eigenvalues have opposite signs, this equilibrium point is a saddle 
point. 

At the equilibrium point S with coordinates (/O3 , -ff|) (see (28)-(29) again), the eigenval- 
ues only change their signs. Thus, point (pg , H^) is also a saddle point — the eigenvalues 
again have opposite signs. 

Finally, at the equilibrium point Q with coordinates {p*2 = Pmin , HI = 0), the eigenvalues 
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are purely imaginary — the characteristic equation is: 

This corresponds to a centre. The trajectories are closed curves, corresponding to a 
cyclic Universe. To analyze what happens and how the temperature parameter affects 
these curves, it would be helpful to expand the dynamical variables p and H near the 
equilibrium point = Pmin -.Hl = 0), namely: p = Pmin + "i" and H = Q + h . The 
dynamical system (21)-(22) becomes: 

^ /i + 3(r - l)/ir + 3aTz(r)/ir2, (40) 



az{T) 



h = -^h^ + ^Tr+^aTz{T)r'^. (41) 
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Figure 3: In linear approximation, the diagram on the phase plane is an ellipse around the centre Q. 
Points R and S are saddles, the origin P attracts trajectories. This is a model of a cyclic Universe. 

Very close to the equilibrium point {p2 = Pmin ■, = 0), i.e. when r and h are very small, 
the linear approximation r = —2>[az{T)]~^h, h = (T/2)r can be easily integrated to 
reveal that r and h sweep an ellipse in the phase plane: 

r = Ccosut, (42) 
h = CQsinut, (43) 
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where C is a constant depending on the initial conditions, 




2az{T) 



3T 



(44) 



is the angular frequency of the oscillations and 




(45) 



is the ratio of the semi-axes of the ellipse. 

It worth noticing that for the cooling Universe, with the drop of the temperature 
(T 0), the angular frequency uj decreases to zero and in result the period of os- 
cillations increases indefinitely and thus periodicity is lost, i.e. only the hot Universe is 
cyclic. At the same time, the ratio Q increases. 

Unlike the case of the origin, where there was a doubly degenerated eigenvalue and the 
degeneracy was dealt with the introduction of the infinitesimal perturbation param- 
eter S, for point (p2 = Pmin , -f^i =0) there is no degeneracy of the spectrum of the 
stability matrix. The stability can be analyzed as follows. The system could be per- 
turbed by variation of the parameters of the model. The purely imaginary eigenvalue 
A, determined in (39), does not pick up, under such parameter perturbation, any real 
part for the following reason. At point = Pmin = 0), the trace of the stability 
matrix L{p^i^ ,0), as seen in (34), is proportional to Hf which does not depend on any 
parameters and is zero. Thus any parameter variation will not alter the trace of the 
matrix L^p^^^ , 0) and in result, the real part of A will always remain zero, as it is half 
of the trace of L{p^i^ ,0). Therefore, the closed trajectories are preserved and stable, 
depending on the initial conditions only. 

It is also very interesting to study the passage of the trajectories on the phase plane 
through regions on the phase plane characterized by inflation. By definition, inflation 
is equivalent to d{t) > and a{t) > 0. Since 



0. Upon substitution of the dynamical equation (11), the condition for inflation becomes 
< —p (the pressure p must therefore be negative). Namely, the regions in the upper 
half [where H > 0, i.e. a{t) > 0] of the phase plane for which H'^ < aTz{T)p'^ — Tp are 
inflationary. Or: 



a • n 

-=H + H^ 

a 



(46) 



and since the scale factor a{t) is always strictly positive, inflation occurs when H+H^ 



> 




(47) 



The boundary of this region is a hyperbola with asymptotes 




(48) 
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and cutting the p-axis exactly at p = Pmin = [0-2(7")]^^ (the left branch of the hyperbola 
is for negative values of p and thus it is of no interest) — see Figure 4. With the drop 
of the temperature (T — )■ 0), the angle CAK between the asymptotes decreases to zero 
and the inflationary region disappears from the phase plane. 



3 

2' 



ff 



-3 

2 4 6 8 10 

P 

Figure 4: The inflationary region is in the upper half-plane to the right of the hyperbola. 

On Figure 4, the equilibrium points are as follows: the origin, the centre F with coordi- 
nates = Pmin , HI = 0), the saddles D and E with coordinates (^3 , H2 3) respectively. 
When p = p% = Pmin(l + 1/^)7 the asymptotes intersect the line DE at points C and 
K with //-coordinates given by zty^oTz^T) (1/2 -|- 1/T) respectively. The saddle 
points D and E are always (i.e. for all T) between the asymptotes — in the inflationary 
region (for D) and symmetrically below it (for E). 

The cyclic Universe (trajectory 1 on Figure 4) experiences an inflationary period when 
the trajectory on the phase-plane is in the upper half-plane to the right of the hyperbola. 
Of course, there is initial data that evolves into open trajectories on the phase plane that 
do not experience inflation at all (trajectories 4, 6, 7), exit a region of initial inflation 
(trajectory 2), enter seemingly eternal inflation (trajectory 3) or stay in a seemingly 
eternal inflation (trajectory 5). In the latter two cases however, with the drop of the 
temperature and the "closure" of the inflationary region, inflation will eventually end. 
Obviously, as the hyperbola goes through the centre F, there is no initial data leading 
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to inflation-only cyclic solutions. 

To sec how much inflation can be "generated" in this model, it is useful to consider 
the e-fold number = Ina. As H = d/a = {d/dt)lna, the e-fold number N is 
equal to J^f H{t)dt. Using the linear approximation (43), see Figure 3, results in 

N = CQsuiuotdt = [CQ /oj) {cos uti — coscutf). During the inflationary phase, 
Luti = and ujtj is approximated by tt/2. Thus the term in the last brackets is approx- 
imately 1. The upper limit of CQ is the i7-coordinate of point D on Figure 4, that 
is = .JpI/S = ^/p~JTTT/T)/^ = ^/{^TT/TY/\3^^^^(T)]. Thus, the upper hmit of 

the e-fold number N is (\/2/3)-\/l/r -f- l/T^, where T is the temperature at the end of 
the inflationary phase. 



4 Van der Waals Gas Revisited 

Many aspects of the van der Waals quintessence scenario have been investigated in a 
series of works of Capozziello et al. [11], [12], [13]. The purpose of this section is to 
briefly recall the properties of the van der Waals quintessece model from the point of 
view of dynamical systems and to compare to the case of virial real gas. 
An interesting property of this model is the fact that the absolute temperature is again 
treated as a parameter, which can be negative. We will retain this peculiarity only in 
view of the fact that a cyclic Universe possibility is realized when the absolute temper- 
ature varies over an entirely negative interval. 

The density of the Universe, p, will be treated as a strictly non-negative quantity. 

In the dynamical equations, the energy density p and the pressure p are that of the van 

der Waals gas. The van der Waals equation of state is [16] (see also [11], [12] and [13]): 

p=Tzrpp-^p^ (49) 

where a = 3pc/pc and /3 = l/{3pc), with pc and pc being the density and pressure of 
the van der Waals gas at the critical point. Here 7 is the temperature of the van der 
Waals gas. It will be treated as a parameter of the model [11], [12], [13]. 
In terms of a dimensionless energy density rj, defined via ij = p/pc, the van der Waals 
equation of state (49) becomes: 

p = - 3pcr?' . (50) 

3 

At the critical point, p = pc (thus rj = 1) and p = pc- At this point, equation (50) 
yileds Pc = {3/8)^pc- Therefore, a = (27/8)7/3 — the three parameters a, (3, and 7 
are not independent [12] and only 7 and pc will be treated as independent parameters. 
Substituting = {3/8)jPc into (50) yields [12]: 
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For the parameter pc the following holds [12]: 

J^c = — , (52) 

Pcrit 

where the quintessence density parameter VLc is equal to (1 — $75 0)/% (here 779 is the 
present day energy density of the van der Waals gas and r^b^o is the baryonic density 
parameter, given by p},{z = 0)/pcrit with z being the red-shift) and the critical density 
of the Universe is pcrit = 3H^/{8ttG) [12]. 

Substituting the van der Waals pressure into the dynamical equations (11) and (12) 
results in: 

fi = -3Hn{l + :^-ljv) = gi{v,H), (53) 

o — T] o 

H = -Ih'-^+3^W^92{v,H), (54) 
2 6 — T] 

where ^ = (3/16)pc = const > is a parameter of the model (instead of the critical 
density pc). 

Next, the existence of equilibrium points will be addressed. Requesting gi{r),H) = 0, 
(53) yields the quadratic equation: 

9 27 

g7r?' - (y 7 + l)r/ + 3(1 + 7) = . (55) 
For real roots to exist, the discriminant 

must be positive. Thus, the range of values of 7, for which two real equilibrium points 
for 77 exist is: 

8 32 ^ 8 32 

-5-45^<^<-5 + 45^ (57) 

or -3.3419 < 7 < 0.1419. 

Again, the physically meaningful range of 7 is: < 7 < — | + ||\/6- 
The corresponding two real values of rj are: 

27 , 1 _t /_135 2 _ 27 , ^ 

Next, 52 (^7, H) = results in four equilibrium values for H: 



^^1*2,3,4(7) 



I2^7jl 



7(3r?r'2 - 9r?r 2 + . 



(59) 



V 3 \ ril,-3 
For these to be real, the following must hold: 

7«2-3)(3r/t;2-H2 + 8)>0. (60) 

As 3?7*^2 - 9^1,2 + 8 = {ril^2 - 3/2)^ + 5/12 > 0, taking into consideration (57), there are 
two regimes: 
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(a) 77^(7), -3.3419 < 7 < 0.1419 (b) 773(7), -3.3419 < 7 < 



(c) »72(7),0 < 7 < 0.1419 



Figure 5: The bigger root 771 (7) [on grapli (a)] is a continuous function of 7, satisfying ?7i (7) < 3 when 
—3.3419 < 7 < and 771(7) > 3 when < 7 < 0.1419. The smaller root 7/2(7) is not continuous [graphs 
(b) and (c)]. It also satisfies 77^(7) < 3 when -3.3419 < 7 < and 773(7) > 3 when < 7 < 0.1419. 
However, in the range — 1 < 7 < 0, 772(7) becomes negative [see graph (b)] and, thus, unphysical. For 
this range of 7, there is only one equihbrium point: 7/1(7). 

(a) -f - i\/6 < 7 < and 7]l 2 < 3 
or 

(b) < 7 < -| + §^/6 and 77^ 3 > 3 

Indeed, when — | — f|\/6 < 7 < 0, 77^2 is always less than 3 and when < 7 < 
— I + ||\/6, r]l2 is always greater than 3: However, for — | — ||\/6 < 7 < —1, there 
is only one positive root: 77^(7) (the bigger root). The other one corresponds to a non- 
physical negative energy density. 

Therefore, there are three intervals for 7, for which equilibrium solutions exist: 

(1) -| - §^/6 < 7 < -1 (i.e. -3.3419 < 7 < -1) 

There are two equilibrium values of r], namely 775}' 2(7), and four corresponding equilib- 
rium values of H, given by (59). 

(2) -1< 7 < 

There is only one positive equilibrium point T/f (7). The corresponding equilibrium values 
of H are now two. 

(3) < 7 < -| + §V6 (i.e. < 7 < 0.1419) 

Again, there are two equilibrium values of rj and four equilibrium values of H. 

At the equilibrium points r]l2{'y) and -^1234(7)) determined from (58) and (59), after 

setting gi(T],H) = and g2{r],H) = in (53) and (54) respectively, the eigenvalues of 



15 



the matrix L{i], H) of the hnearized dynamical system are: 

A. = (61) 

A2 = -3ffi*2.3,,(7). (62) 

As the physicahy relevant 77* 2(7) are positive [r] is energy density), the signs of Ai^2 are 
determined by the sign of 7, the sign of -^1234(7)) ^-i^d the sign of the expression in 
the figure brackets in (61). 

The latter term, [3 - ryj 2(7)]"^ - 3/8, is positive when: 3r/j|" 2(7) - 18r/j|" 2(7) + 19 < or 
3 - (2/3)V^ < r?|,2(7) < 3 + (2/3)^/6. That is, 1.3670 < 77^,2(7) < 4.6330. 
In the first interval of 7 for which equilibrium points exist, 

(1) -| - § V6 < 7 < -1 (i.e. -3.3419 < 7 < -1), 

there are four equilibrium points: 
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Fi gure 6: For — 3.3419 < 7 < — 1 (and ^ = 1), points A and B are saddles, point C — stable node, 
point D — unstable node. 



(i) Point A with coordinates -^1(7) > and 1.3670 < 7/^(7) < 4.6330. At this point 
Ai > and A2 < 0. This is a saddle point. 
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(ii) Point B with coordinates H^i^f) < and 1.3670 < r]l{'y) < 4.6330. The eigenval- 
ues flip signs: Ai < and A2 > 0. Another saddle point. 

(iii) Point C with coordinates -^3(7) > and ??2(7) < 1-3670. At this point, both 
eigenvalues Ai_2 are negative and this is a stable node. 

(iv) Point D with coordinates HKj) < and r/2(7) < 1-3670. At this point, both 
eigenvalues Ai_2 are positive, corresponding to an unstable node. 

In the next interval of 7 with equilibrium points: 
(2) -1< 7 < 0, 

there are only two equilibrium points (recall that ?/2(7) negative in that range): 
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Figure 7: For — 1 < 7 < (and ^ = 1), there are only two equilibrium points — E and F — both 
saddles. 



(i) Point E with coordinates i^i(7) > and 1.3670 < rjl{-i) < 4.6330. At this point 
Ai > and A2 < — saddle point. 

(ii) Point F with coordinates i^|(7) < and 1.3670 < 7?J(7) < 4.6330. The eigenval- 
ues flip signs: Ai < and A2 > — saddle point. 
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This is again a cyclic Universe (but corresponding to a toy model with negative absolute 
temperature). 

Finally, in the last interval of 7 with equilibrium points: 

(3) < 7 < -| + §VG (i.e. < 7 < 0.1419), 
there are again four equilibrium points: 
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Figure 8: For < 7 < 0.1419 (and ^ = 1), points M and A*' are saddles, point G — stable node, 
point K — unstable node. 

(i) Point G with coordinates H^{j) > and 1.3670 < 77^(7) < 4.6330. Both eigen- 
values Ai^2 are negative and this is a stable node. 

(ii) Point K with coordinates -ff|(7) < and 1.3670 < 772(7) < 4.6330. Both eigen- 
values Ai^2 are positive and this is an unstable node. 

(iii) Point M with coordinates -^^3(7) > and 4.6330 < 7/^(7). At this point Ai > 
and A2 < 0. This is a saddle point. 

(iv) Point N with coordinates -^^4(7) < and 4.6330 < rjK'j). The eigenvalues flip 
signs: Ai < and A2 > 0. Again, a saddle point. 

The point (r/* = 0, H* = 0) in the phase portrait is also an equilibrium point. Appar- 
ently, at this point, Ai = = A2 is a double eigenvalue of the stability matrix L(0, 0) — 
as in the case of the virial real gas. Again, the introduction of an infinitesimal parameter 



18 



5 > in the same way in H perturbs the dynamical system and bifurcates the origin 
into a "dipole" of two equihbrium points. We consider: 

7? = -3il77 (1 + - ^7^) , (63) 
6 — 1] 8 

H = -^H' + ls'-^^+S^W. (64) 

Clearly, {rj^ = 0,Hq* = +6) and (773 = 0,Hq* = —5) are the two equilibrium points of 
the "dipole". At these points, the eigenvalues of the stability matrix L(0, ±5) are: 

Ai = -3(±5)(l + 7), (65) 
A2 = -3{±5). (66) 

There arc only two possible cases — the situation where 7 is positive (which leads to 
equilibrium points with ry* greater than 3) cannot be studied at the origin. 
Firstly, if — 1 < 7 < 0, then Ai and A2 have the same signs. Thus, point (77^ = 0, H^* = 
+6) is a stable node, while point {rj^ = 0, H^* = —S) is an unstable node. In the limit 
(5—7-0, the two points coalesce and this reveals that the origin attracts trajectories in 
the phase plane. The graph on Figure 7 is realized. 

The other possible case is —3.3419 < 7 < —1. Then Ai and A2 have the different signs 
and the two equilibrium points of the "dipole" are saddle points. This shows that, upon 
taking the limit (5 — >■ 0, the origin will not be able attract any trajectories in the phase 
plane. 

5 Conclusions 

It has been shown that there is a possible cyclic Universe scenario with inflation real- 
izable within a classical physics model which offers a graceful exit from inflation. The 
cyclic solution is possible for a wide range of initial data and a wide range of tem- 
peratures for real gas with equation of state derived from the virial expansion and 
for negative temperatures for a van der Waals real gas. Periodicity is lost too as the 
temperature of the Universe tends to zero. 
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